Hybrid method for understanding black-hole mergers: Inspiralling case 
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We adapt a method of matching post-Newtonian and black-hole-perturbation theories on a time- 
like surface (which proved useful for understanding head-on black-hole-binary collisions) to treat 
equal-mass, inspiralling black-hole binaries. We first introduce a radiation-reaction potential into 
this method, and we show that it leads to a self-consistent set of equations that describe the simul- 
taneous evolution of the waveform and of the timelike matching surface. This allows us to produce 
a full inspiral-merger-ringdown waveform of the I = 2, m = ±2 modes of the gravitational wave- 
form of an equal-mass black-hole-binary inspiral. These modes match those of numerical-relativity 
simulations well in phase, though less well in amplitude for the inspiral. As a second application 
of this method, we study a merger of black holes with spins antialigned in the orbital plane (the 
superkick configuration). During the ringdown of the superkick, the phases of the mass- and current- 
quadrupole radiation become locked together, because they evolve at the same quasinormal-mode 
frequencies. We argue that this locking begins during the merger, and we show that if the spins of the 
black holes evolve via geodetic precession in the perturbed black-hole spacetime of our model, then 
the spins precess at the orbital frequency during the merger. In turn, this gives rise to the correct 
behavior of the radiation, and produces a kick similar to that observed in numerical simulations. 

PACS numbers: 04.25.Nx, 04.30.-w, 04.70.-s 



I. INTRODUCTION 

Black-hole-binary mergers are both key sources of 
gravitational waves [1| and two-body systems in general 
relativity of considerable theoretical interest. It is com- 
mon to describe the dynamics and the waveform of a 
quasicircular black-hole binary as passing through three 
different stages: inspiral, merger, and ringdown (see, e.g., 
0). For comparable- mass black holes, the three stages 
correspond to the times one can use different approxi- 
mation schemes. During the first stage, inspiral, the two 
black holes can be modeled by the post-Newtonian (PN) 
approximation as two point particles (see, e.g., Q for a 
review of PN theory) . As the speeds of the two holes in- 
crease while their separation shrinks, the PN expansion 
becomes less accurate (particularly as the two objects 
begin to merge to form a single body). In this stage, 
merger, gravity becomes strongly nonlinear (and there- 
fore less accessible to approximation techniques). Af- 
ter the merger, there is the ringdown, during which the 
spacetime closely resembles a stationary black hole with 
small perturbations [and one can treat the problem using 
black- hole perturbation (BHP) theory (see, e.g., [1] for a 
review of BHP theory)]. 

Because the merger phase of comparable-mass black 
holes has been so challenging to understand analytically, 
there have been many attempts to study it with a vari- 
ety of analytical tools. One approach has been to develop 
PN and BHP theories to high orders in the different ap- 
proximations. Since neither approximation can yet de- 
scribe the complete merger of black-hole binaries, several 
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groups worked on developing methods that aim to get 
the most out of a given approximation technique. The 
close-limit approximation (see, e.g., for early work 

and for more recent work) and the Lazarus project 

(see, e.g., [ll,|i3) both try to push the validity of BHP 
to early times; the effective-one-body (BOB) approach 
(see, e.g., [H, [l^ for the formative work, and |17l - [20| 
for further developments that allow the method to repli- 
cate numerical-relativity waveforms) aims to extend the 
validity of the PN approximation to later times. 

There also have been several methods that do not 
easily fit into the characterization of extensions of 
PN or BHP theories. For example, the "particle- 
membrane" approach of Anninos et al. [2ll, [24I com- 
putes the waveform from head-on collisions by ex- 
trapolating results from the point-particle limit to the 
comparable-mass case (and taking into account changes 
to the horizons computed within the membrane paradigm 
^23*1). More recently, white- hole fission was used in ap- 
proximate models of black-hole mergers |2J- T26:j . and 
quite recently, Jaramillo and collaborators |27H30| used 
Robinson- Trautman spacetimes as an approximate ana- 
lytical model of binary mergers (as part of a larger project 
correlating geometrical quantities on black-hole horizons 
with similar quantities at future null infinity). 

Analytical approximations are not limited to 
comparable-mass black-hole binaries, and recently 
there has been a large body of work on developing 
techniques to study intermediate- and extreme-mass- 
ratio inspirals (IMRIs and EMRIs, respectively). Most 
of these methods aim to produce gravitational waves 
in ways that are less computationally expensive than 
computing the exact numerical solution or computing 
the leading-order gravitational self- force are (see, e.g., 
[sH for a recent review of the self- force). The majority 
of the approaches rely heavily on BHP techniques 
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combined with some prescription for taking radia- 
tive effects into account, though not all approximate 
methods fall into this classification (Barack and Cutler 
[3^ . for example, model EMRIs by instantaneously 
Newtonian orbits whose orbital parameters vary slowly 
over the orbital time scale because of higher-order PN 
effects). A well-known example is that of Hughes [33| . 
Glampedakis 34], Drasco ^35,], Sundararajan [3^ and 
their collaborators whose semi-analytical approaches are 
often called Teukolsky-based models. These methods 
describe the small black hole as moving along a sequence 
of geodesies whose energy, angular momentum, and 
Carter constant change from the influence of emitted 
gravitational waves. They usually involve some ad- 
ditional prescription to treat the transition from the 
inspiral to the plunge, when the motion is no longer 
adiabatic. The EOB formalism in the EMRI limit, 
however, does not require an assumption of adiabatic 
motion (see, e.g., [37-42]). By choosing the dynamics 
of the EMRI to follow the EOB Hamiltonian and a 
resummed multipolar PN radiation- reaction force [43| . 
these authors can calculate an approximate waveform 
without any assumption on relative time scales of orbital 
and radiative effects. One can also make an adiabatic 
approximation with EOB methods, as Yunes et al. 
(44 . I45I recently did in their calibration of the EOB 
method to a set of Teukolsky-based waveforms. Lousto 
and collaborators (46l - l48| took a different approach to 
the EMRI problem in their recent work. They used tra- 
jectories from numerical-relativity simulations of IMRIs 
as a way to calibrate PN expressions for the motion of 
the small black hole. They then performed approximate 
calculations of the gravitational waves using the PN 
trajectories in a black-hole perturbation calculation, and 
found good agreement with their numerical results. 

In a previous article [i^ (hereafter referred to as Pa- 
per I), we showed that for head-on collisions, one can 
match PN and BHP theories on a timelike world tube 
that passes through the centers of the PN theory's point 
particles. The positions of the points particles as a func- 
tion of time (and, consequently, the world tube) were 
chosen before evolving the waveform. Moreover, they 
were selected in such a way that both PN and BHP the- 
ories were sufficiently accurate descriptions of the space- 
time on the world tube or the errors in the theories did 
not enter into the waveform. (A plunging geodesic in the 
Schwarzschild spacetime worked in Paper I.) This allowed 
us compute a complete waveform for all three phases of 
black-hole-binary coalescence and gave us a way to in- 
terpret the different portions of the waveform. More- 
over, when we compared the waveform from the hybrid 
method with that of a full numerical simulation of plung- 
ing equal-mass black holes with transverse, antialigned 
spins, we found very good agreement between the two. 

There is no reason, a priori, why the same procedure 
of Paper I (namely, specifying the position of the point 
particles as a function of time and matching the metrics 
on a surface passing through their positions) should not 



work for inspiralling black holes as well. The principal 
difficulty arises from trying to find a way of specifying the 
positions of the particles for inspiralling black holes (and 
thus a location at which to match the PN and BHP met- 
rics) that does not introduce errors into any of the three 
stages of the inspiral, merger, or ringdown portions of the 
waveform. The most important development that we in- 
troduce in this paper, therefore, is a way of achieving this 
goal by including a radiation-reaction force into the for- 
malism. In the hybrid method, we compute a radiation- 
reaction force by using the outgoing waves in the exterior 
BHP spacetime to modify the PN dynamics in the inte- 
rior through a radiation-reaction potential ^] . We show, 
in this formalism, that introducing a radiation-reaction 
potential is equivalent to solving a self-consistent set of 
coupled equations that describe the evolution of the point 
particles' reduced-mass motion and the outgoing gravita- 
tional radiation, where the particles generate the metric 
perturbations of the gravitational waves and the waves 
carry away energy and angular momentum from the par- 
ticles (thereby changing their motion). 

Our principal goal in the paper is to explore this cou- 
pled set of evolution equations and show, numerically, 
that it gives rise to convergent and reasonable results. 
We will use these results to make a refinement of our 
interpretation of the waveform from Paper I, and we 
will also compare the waveform generated by the hybrid 
method to that from a numerical-relativity simulation of 
an equal-mass, nonspinning inspiral of black holes. The 
two waveforms agree well during the inspiral phase, but 
less well during merger and ringdown. The discrepancy 
at late times is well understood: we continue to model 
the final black hole produced from the merger as nonspin- 
ning, although, in fact, numerical simulations have shown 
the final hole to be spinning relatively rapidly (see, e.g., 
[5l|). Adapting the hybrid approach to treat the final 
black hole as rotating is beyond the scope of this work, 
but is something that we will investigate in the future. 

As an application of the hybrid method for inspirals, 
we explore the large kicks produced from black-hole bi- 
naries with antialigned spins in the orbital plane (the su- 
perkick configuration [52|, |53|). As noted by Schnittman 
et al. [S^l and emphasized to us by Thorne , the spins 
must precess at the orbital frequency during the final 
stage of the merger. While Briigmann et al. !56] were 
able to replicate this effect using a combination of PN and 
numerical-relativity results, we will need to take a differ- 
ent approach, by using geodetic precession in the exterior 
Schwarzschild BHP spacetime, to have the spins lock to 
the orbital motion at the merger. When we include the 
geodetic effect, we are able to recover the correct qual- 
itative profile of the kick, although the magnitude does 
not match precisely. 

We organize the paper as follows: We review the re- 
sults of Paper I in Sec. [Hi and we describe the procedure 
for calculating the radiation-reaction force and the re- 
sulting set of evolution equations in Sec. IIIIl In Sec. IIVI 
we show the convergence of our waveform, we compare 
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Matching Shell BMP r=R+M 

FIG. 1: (Reproduced from Paper I.) The regions of spacetime 
and the radial coordinates in the hybrid method (at a given 
moment in time, with one spatial coordinate suppressed). The 
exterior is a perturbed Schwarzschild and the interior is a PN 
spacetime. At the position of the shell, the two descriptions of 
spacetime should both be valid, or should be within a region of 
spacetime that does not heavily influence physical observables 
far away. 



with numerical relativity, and we discuss using the hybrid 
method to interpret the waveform. Next, we discuss the 
behavior of spinning black holes and describe spin pre- 
cession as a mechanism for generating large black-hole 
kicks in Sec. El We conclude in Sec.|VTl Throughout this 
paper, we set G — c — 1, and we use the Einstein sum- 
mation convention (unless otherwise noted). 



II. A BRIEF REVIEW OF PAPER I 

In this section, we will review the essentials of the for- 
malism from Paper I. In the hybrid method, we divide 
the spacetime of an equal-mass, black-hole-binary merger 
into two regions: a PN region within a spherical shell 
through the centers of the PN theory's point particles, 
and a perturbed Schwarzschild spacetime outside that 
shell. Figure [1] shows this at a given moment in time 
(with one spatial dimension suppressed). For the hybrid 
procedure to work, there must be either a spherical shell 
on which both BHP and PN theories are simultaneously 
valid (to a given level of accuracy) or a way to prevent 
the errors in the approximations from affecting observ- 
ables, such as the waveform. By finding good agreement 
between the hybrid waveform and that of numerical rel- 
ativity in Paper I, we found evidence that matching the 
theories on a spherical shell that passes through the PN 
theory's point particles works throughout all three stages 
of a head-on black-hole-binary merger: infall, merger, 
and ringdown. 

To mesh the two descriptions of spacetime, we match 
the PN metric to that of the perturbed Schwarzschild 
black hole, which involves relating the two coordinate 
systems of PN and BHP theories. In the PN coordi- 
nate system, we will use uppercase variables, and we 
will use a harmonic gauge. For example, we will employ 
(T, X, Y, Z) when describing the Cartesian coordinates of 
the background Minkowski space and {T,R, 8, $) when 
discussing its spherical-polar coordinates. In the per- 
turbed Schwarzschild spacetime, we will use (t,r,9,ip), 
primarily, though sometimes we will also use the light- 



cone coordinates, {u,v,9,ip), where 

u = t — r^, , V = t + r^, , 

and 



2M log 



2M 



- 1 



(1) 



(2) 



One can match the two coordinate systems, accurate to 
linear order in M / R by identifying 

T = i, 9 = 61, $ = R = r-M . (3) 

For the equal-mass binaries that we study, we will de- 
note the separation by A{t) = 2R(t) in PN coordinates 
and a{t) = 2r{t) in Schwarzschild coordinates. More- 
over, because we match the two metrics on a shell pass- 
ing through the centers of the point particles, we will 
indicate the position of the shell by adding a subscript 
"s" to the coordinate radius. For example, we will write 
R^{t) = A{t)/2 or rs{t) = a{t)/2 to denote this. For 
clarity, we reproduce a table that reviews the essentials 
of our notation in Table H] 

Because we are investigating only the lowest-order ef- 
fects in our study of radiation reaction and large black- 
hole kicks, we shall only need the lowest-order terms in 
the PN metric that appeared in Paper I to describe the 
interior of the shell. 



dS^ = -(1 - 2M/i?- 2C/^-^Vt2 

(;=2)s 



(;=2) 



dtdx" 



(1 + 2M/R + 2c/;;-^0(rf^ + R d'n) 



(4) 



In the above equation, M is the total mass of the binary, 
d^H. is the area element on the unit sphere, dx^ = d9, dip, 
and the additional variables J/]^"^'' and w^'"^"* are the 
quadrupole parts of the spherical harmonic expansion 
of the binary's Newtonian potential and gravitomagnetic 
potential, respectively. 



= E u'/''Y,,^{e,^), (5) 

m=-2 
2 



W 



(1=2) 



(6) 



m=-2 



We denote the scalar spherical harmonics by l2.m(6', </?), 
and the coefficients J/^™ and w^^™ are functions of R 

and t. The functions X^'™(0,(^) are odd-parity vector 
spherical harmonics, whose 9 and ip coefficients are given 
by 



= -{csc9)d^Y'^-"\9,^) 



(7) 

(8) 



A more general description is put forth in Paper I, but 
here we only take the essential components needed for 
the calculations in the paper. 
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PN spacetime 


Matching 


; shell 


Perturbed Schwarzschild spacetime 


Coordinates 




it,Rs{t),e,if) or 


{t,r,{t),e,ip) 


{t,r,e,<fi), r = R + M 


Binary separation 


A{t) 


A{t) or 


a{t) 


ait) 


Matching radius 


R{t) = A{t)/2 


Rs{t) = a{t)/2 - M or 


rs{t) = A{t)/2 + M 


r{t) = a{t)/2 



TABLE I: (Reproduced from Paper I.) The notation for the coordinates, the binary separation, and the matching radius. We 
express these three variables in the PN spacetime, the BHP spacetime, and the matching surface between the two. 



Outside of the shell, we write down a perturbed 
Schwarzschild metric, 

ds^ = -(1 - 2M/r)dt^ + (1 - 2M/r)-\dr^ + r^d'^n) 
+h^^dx^'dx'', (9) 

where the nonzero components of the perturbed metric 
/i^ii, that we shall need in this paper are the quadrupole 

pieces, /i^'i7^'' , and they take the form, 

(ftir^e) = E Hfrr^'-^ie,^), (10) 

m=-2 
2 

e) = E H^rY^'^'iO,^), (11) 

(12) 



m=-2 
2 



m=-2 



sm^e J2 K^^"'Y^^"\0,^), (13) 



m=-2 



and 



- E h'rxf"\o,^), (14) 



(Kv )(o) 



m=-2 
2 



E /.?'™x2-(^,^). (15) 



The subscripts (e) and (o) refer to the parity of the 
perturbations (even and odd, respectively), where we 
call perturbations that transform as (—1)' even and as 
(-1)'+! odd. 

The interior PN metric must match the perturbed 
Schwarzschild metric on a spherical shell between the 
two regions. To make this identification, we note that 
because R~r — M, then the term 



+ j + 0[{M/Rf] . (16) 



We, therefore, identify the monopole piece of the PN met- 
ric with the unperturbed Schwarzschild metric. More- 
over, at leading order in M / R, we note that the pertur- 
bations of the two metrics match exactly. 



2,m 



N ' 



-Aw 



(o) 



(17) 
(18) 



There is then a straightforward procedure that lets one 
express the metric perturbations in terms of the gauge- 
invariant perturbation functions of the Schwarzschild 
spacetime [svj (though in this paper we use the nota- 
tion of [H^), which are typically called the Zerilli func- 
tion and the Regge- Wheeler function for the even- and 
odd-parity perturbations, respectively. We reproduce the 
expressions below: 



2,m 



2r f 2,m , r~2M 
3 ^ 2r + 3M 



1 - 



2M 



U 



2,m 
TV 



rdrUj^ 



2,m 



(o) 



2r ( drW 



(o) 



^ 2,m 
-W,\ 

r (°) 



(19) 



(20) 



In Paper I, we matched the two metrics on a timelike 
tube that we specified before evolving the Regge- Wheeler 
and Zerilli functions. We assumed that this tube would 
be spherically symmetric, and we found its radius by first 
assuming the reduced-mass motion of the system followed 
a radial geodesic of a plunging test mass in the back- 
ground Schwarzschild spacetime and then setting the ra- 
dius of the world tube to be half this distance at each 
time. This allowed us to use the PN data in the form of 
the Regge- Wheeler and Zerilli functions, Eqs. ((TO)) and 
(l20l) on this tube to provide a boundary-value problem 
for the evolution of the Regge- Wheeler [s^ or Zerilli [60l | 
equations. 



(c,o) 

dudv 



V,o)'^(o,o) 



0. 



(21) 



The potentials for the Regge- Wheeler (odd-parity) or 
Zerilli (even-parity) equations are given by 



^(c,o)(0 



2M 



^C/' (r) 

Tr^(o,o)W 



(22) 



where A = ^(/ + 1) and 

C/(',)(r) = l, (r) 



A(A -I- 2)r2 -f 3M(r - M) 



(Ar + 3M)2 

(23) 

where A ^ (l - + 2)/2 = A/2 - 1. After numerically 
solving the Regge- Wheeler or Zerilli equations above, we 
computed the gravitational waveforms and the radiated 
energy and momentum, all of which we found to be in 
good agreement with the exact quantities computed from 
numerical- relativity simulations. 
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In this paper, while much of the procedure we use for 
matching the metrics is identical to that set forth above, 
there are several important differences that we will dis- 
cuss in Sec. IIIII The most important difference between 
the first paper and the current one arises in how we find 
the trajectory of the system's reduced mass (and then the 
timelike tube on which we match the metrics). Before, 
we chose a region, prior to evolving the Regge- Wheeler 
and Zerilli equations, that would not introduce spurious 
effects into the results; here we determine the position 
of timelike world tube through evolving the position of 
the reduced mass of the binary subject to a radiation- 
reaction force. We will discuss the details of this proce- 
dure in the next section. 



III. RADIATION-REACTION POTENTIAL 
AND EVOLUTION EQUATIONS 

In this section, we introduce a radiation-reaction po- 
tential into the hybrid method, and we show that it leads 
to a set of evolution equations that simultaneously evolve 
both the outgoing radiation and dynamics of the reduced 
mass of the system. This, in turn, allows us to produce 
a full inspiral-merger-ringdown waveform. We first qual- 
itatively discuss how our method works and how it com- 
pares to other analytical methods. We then discuss the 
hybrid method in further detail, and we close this section 
by showing, analytically, that the procedure recovers the 
correct Burke-Thorne radiation-reaction potential f50| in 
the weak- field limit. 



A. Qualitative Description 

It is easiest to discuss our method with the aid of 
the spacetime diagram in Fig. [21 We describe the re- 
gion within the solid black timelike curve with the near- 
zone PN metric, and outside this curve, we use a per- 
turbed Schwarzschild region. The black line, which 
passes through the PN point particles, is where we match 
the two metrics. We suppress both angular coordinates, 
so that each point on the curve represents the matching 
shell that we discuss in Sec. HIl The shaded region rep- 
resents the black-hole potential; the yellow (light gray) 
shade depicts the strong-field portion of the black-hole 
potential (the strong-field near zone) and the green (gray) 
shade shows the region where centrifugal potential is sig- 
nificant (the weak- field near zone) . There is a wave zone 
near the horizon (large u), and, consequently, the region 
where there is a large black-hole potential is confined to 
a small space in this diagram. 

To have Fig. [5] be an effective description of the space- 
time of a black-hole binary, both PN and BHP theories 
both must be sufficiently accurate at the PN point parti- 
cles (the black line where we match the metrics), or the 
PN approximation could break down if the point particles 
are well-hidden within the black-hole effective potential. 




FIG. 2: A spacetime diagram of our method. The solid black 
timehke curve depicts the region where we match PN and 
BHP spacetimes (passing through the centers of the PN the- 
ory's point particles). Inside this cmve, the spacetime can 
be reasonably approximated by PN theory, whereas outside, 
the spacetime is better described by BHP theory. The yel- 
low (light gray) shade shows the strong-field region, whereas 
the green (gray) shade represents where the black-hole po- 
tential is weaker, but the centrifugal barrier of flat space still 
is important. The dark blue (dark gray) shaded region, sur- 
rounded by the blue (dark gray) dashed lines shows how the 
value of the perturbations in the exterior (along with the no- 
ingoing- wave condition) determines the value of the radiation- 
reaction potential at the next matching point. The horizon- 
tal dashed lines represent the region of spacetime where the 
close-limit approximation or Lazarus approach would begin. 
The red (gray) dashed lines show how one can connect the 
near-zone behavior to the wave (through lines of constant u), 
and thereby tie the motion of the matching region through 
the black-hole effective potential to portions of the waveform. 
This gives an interpretation of inspiral, merger, and ringdown 
phases in terms of the direct and scattered parts of the waves. 
Further discussion of this figure is given in the text of Sec. IIIII 



For the errors to stay within the potential, the particles 
must rapidly fall to the horizon; thus one can see in Fig. 
[5] that the black curve approaches an ingoing null ray 
asymptotically. (Recall that u and v are the light-cone 
coordinates of the BHP spacetime.) Following this tra- 
jectory, the perturbations induced by the PN spacetime 
will become strongly redshifted, and they will not escape 
the black-hole potential (as Price had found in his de- 
scription of stellar collapse |61J), because the potential 
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reflects low frequency perturbations. 

As in Paper I, we will again be able to interpret differ- 
ent portions of the waveform by connecting a region of 
the waveform with the position of the PN binary's point 
particles in the near zone (via constant values of the light- 
cone coordinate, u). In the figure these are the thick red 
(gray) dashed lines of constant u. The inspiral part of the 
waveform, which propagates directly along the light cone, 
comes from the part of the trajectory within the weak- 
field near zone. Once the trajectory reaches within the 
strong-field near zone, the waves scatter off of the poten- 
tial and propagate within the light cone (often referred to 
as the PN tail part of the wave) in addition to propagat- 
ing out directly. We view this mixed wave as character- 
istic of the merger phase. Finally, as the trajectory falls 
within the effective potential, for a Schwarzschild black 
hole the direct part vanishes and only the scattered waves 
emerge; this part is the quasinormal ringing of the final 
black hole and should be associated with the ringdown 
phase. We distinguish between Schwarzschild and Kerr 
black holes for the ringdown phase, because Mino and 
Brink (62i] and, subsequently, Zimmerman and Chen [GSJ 
found that for Kerr black holes, frame dragging generates 
a part of the waveform at the horizon frequency (that de- 
cays at a rate proportional to the horizon's surface grav- 
ity). This piece of the waveform looks like a source, and, 
thus, only when the final black hole is not spinning do we 
consider the spacetime to appear to be source free. For 
this reason (and since the matching surface asymptotes 
to a line of constant v), we call the values of v greater 
than this limiting value the homogeneous region, and the 
values of v less than this the source region. 

An important development in this paper is that we no 
longer prescribe the evolution of the reduced mass of the 
system (and thereby a matching region) before evolving 
the Regge- Wheeler or Zerilli equations; rather, we spec- 
ify a set of evolution equations for the conservative dy- 
namics of the binary, and let the outgoing waves provide 
back reaction onto the dynamics. This, in turn, leads 
to a self-consistent system of equations including radia- 
tion reaction. More concretely, we continue to match the 
PN and perturbed Schwarzschild metrics at the centers 
of the PN theory's point particles. Moreover, we will 
again let the reduced-mass motion of the binary system 
follow that of a point particle in a Schwarzschild back- 
ground; in this paper, however, we will use the fact that 
there are no ingoing waves to specify a radiation-reaction 
potential that acts as a dissipative force on the Hamil- 
tonian dynamics of the reduced mass. This follows the 
spirit of the Burkc-Thorne radiation-reaction potential, 
but the radiation propagates within a BHP spacetime, 
and, therefore, also takes the effects of the background 
curvature into account. 

Furthermore, adding a radiation-reaction force to the 
hybrid method leads to a set of equations that simultane- 
ously evolve the Zerilli equation (the waveform) and the 
reduced-mass motion of the binary. In Fig. [5] we repre- 
sent schematically how this occurs. We start at a given v 



[a dark blue (dark gray) dashed line given in Fig. [2] and 
assume that there is a no-ingoing-wave boundary condi- 
tion along the line u — 0. In addition, we suppose that we 
have determined the black-hole-perturbation functions 
for all smaller values of v, up to the timelike match- 
ing surface. By evolving the Zerilli equation, Eq. (|21l) . 
one can find the Zerilli function at v + dv up to the 
time Us{v) [within the dark blue (dark gray) shaded re- 
gion] . The no-ingoing- wave condition combined with the 
boundary condition on the matching surface, however, 
fixes how the Zerilli function will evolve to larger values 
of u. When solved simultaneously with the Hamiltonian 
dynamics describing the binary's motion, this lets one 
find the position of the reduced mass of the binary at 
V + dv, denoted by Us{v + dv) in the figure, and the new 
value of the Zerilli function there. One can evolve the 
system for all v in such a manner. 

Including a radiation-reaction force does not greatly 
change the hybrid method as reviewed in Sec. [TTl The 
matching procedure works the same; one modification 
that comes about is that we must include both the 
Newtonian potential and the radiation-reaction poten- 
tial in the PN metric (and, therefore, gain an additional 
term in the Zerilli function). The evolution system is 
now quite different, because it is a coupled system of 
Hamiltonian ordinary differential equations and a one- 
dimensional partial differential equation. We will discuss 
the system of evolution equations in greater detail after 
we compare our method with other analytical methods 
in the next subsection. 



B. Descriptive Comparison with Other Analytical 
and Semi-Analytical Models 

In this section, we will compare the similarities and 
differences between the hybrid method described above 
and the most closely related methods mentioned in 
the Introduction: the close-limit approximation, the 
Lazarus program, the comparable-mass EOB meth- 
ods, the Teukolsky-based approach, the EOB descrip- 
tion of EMRIs, and the IMRI calculation calibrated to 
numerical-relativity data. The comparison between the 
hybrid method and the other methods will be descrip- 
tive, but we will compare the waveform from the hybrid 
method with a numerical-relativity waveform in Sec. IIVI 

To compare with the Lazarus project or the close-limit 
approximation, we again refer to Fig. [31 where we show 
two spacelike hypersurfaces (the horizontal dashed lines 
labeled by Ei and ^,2). In the close-limit and Lazarus 
methods, initial data is posed on these surfaces at a time 
near the merger of the black holes. While these ap- 
proaches have been successful, posing initial data at late 
times makes it more difficult to smoothly connect the 
initial inspiral of the binary to the merger and ringdown 
later. Moreover, because the initial data extends inside 
the black-hole potential, if it contains high-frequency per- 
turbations, these could escape the potential barrier and 
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enter into the waveform. The hybrid approach escapes 
this problem by setting boundary data on a timelike 
world tube rather than on a spacelike hypersurface. This 
also lets the method connect the inspiral, merger, and 
ringdown portions of the dynamics and waveform more 
directly. 

The EOB approach, for comparable-mass ratio bina- 
ries, only describes times prior to the merger (the hyper- 
surfaces Si and E2 in Fig. [2|). To create a full inspiral- 
merger-ringdown waveform, the EOB method must fit a 
sequence of quasinormal modes to the end of the insprial- 
plunge waveform. This procedure makes a very accurate 
waveform, but it makes connecting the behavior of the 
spacetime before and after the merger more difficult. The 
hybrid method, with its interior PN region that falls to- 
ward the horizon at late times, allows one to make a more 
clear connection between the dynamics of the spacetime 
during inspiral and merger to that during ringdown. In 
its current implementation, however, it does not produce 
a waveform nearly as accurate as that of the EOB. 

Although the hybrid method is designed for describ- 
ing comparable-mass black-hole binaries, it shares a few 
similarities and has several significant differences from 
various approximate techniques that model EMRIs. It 
is possible to draw a few general comparisons between 
the hybrid method and the procedures for studying EM- 
RIs, before moving to more specific comparisons. While 
the hybrid method evolves perturbations on a black-hole 
background (as most EMRI methods do), EMRI methods 
assume a source term as the generator of the perturba- 
tions in the background. The hybrid approach, however, 
does not have a source term; rather, the perturbations 
of the background come from boundary data that corre- 
spond to the multipolar structure of a comparable-mass 
PN binary. Because the hybrid approach is a boundary- 
value problem, the details of the implementation will be 
different from those methods that use a point mass as a 
source term. 

Moving to specific EMRI models, we first compare 
the hybrid approach with the Teukolsky-based meth- 
ods of Sundararajan and collaborators [36| (for exam- 
ple). The hybrid approach is similar to that of [36| , 
in that both use time domain codes and are capable of 
producing smooth inspiral-merger-ringdown waveforms. 
An important difference is that the hybrid method cal- 
culates the waveform simultaneously with the evolution 
of the matching region, whereas the EMRI method of 
Sundararajan computes the trajectory before the evolu- 
tion (using an adiabatic frequency-domain code during 
insprial, and a prescription for the plunge and merger) 
and then finds the waveform from this trajectory. More- 
over, we compute the radiation-reaction force in the hy- 
brid method by matching the near-zone PN solution to an 
outgoing solution in the exterior BHP spacetime, whereas 
the Teukolsky-based methods include radiative effects by 
evolving the orbital parameters of geodesies from aver- 
aged fluxes at infinity. 

The EOB model of Yunes et al. {i^, is a calibration 



of the EOB method to Teukolsky-based waveforms for 
EMRIs; it, therefore, shares the same similarities and 
differences as the EOB and the Teukolsky-based meth- 
ods discussed above. Han and Cao [i^ develop an EOB 
model that uses the a Teukolsky-based energy flux (in the 
frequency domain) to treat radiative effects. In compar- 
ing with the hybrid model, therefore, it also falls some- 
where between an EOB model and a Teukolsky-based 
method. The recent EOB work of Bernuzzi and collab- 
orators [39I - I4H shares more similarity with the hybrid 
method, because they evolve the Regge- Wheeler- Zerilli 
equations in the time domain. The most notable specific 
difference (as opposed to the general differences between 
the hybrid-method and all analytical approaches to EM- 
RIs noted above) is in the radiation-reaction force. The 
EOB model uses a high-PN-order, resummed energy flux, 
whereas (as also noted above) the hybrid method deter- 
mines radiative effects from directly matching a near- 
zone PN solution to an outgoing BHP solution. 

We conclude this section by comparing the hybrid 
method with the recent analytical work of Lousto et 
al. [46l - l48j . They take two approaches to calculating 
waveforms for IMRIs perturbatively. In their initial 
work, they transform the trajectory of the small black 
hole from their numerical-relativity simulations into the 
Schwarzschild gauge, and they compute the waveform us- 
ing this numerical trajectory in a BHP calculation. To 
be able to study a wider range of mass ratios, they use 
PN expressions for the change in frequency and the ra- 
dial trajectory, but use the numerical-relativity values of 
the frequency to calibrate the PN functions. The hy- 
brid method differs from this, because it calculates the 
matching region simultaneously with the waveform, and 
it does not use numerical-relativity data to calibrate re- 
sults. Consequently, the hybrid method does not agree 
as well with exact results as well as the other methods 
discussed here, but it does present a distinct way of calcu- 
lating the approximate spacetime and gravitational wave- 
form. 



C. Radiation Reaction and Evolution Equations 

In this section, we will discuss the details of radiation 
reaction in the hybrid method. The end result will be 
the set of evolution equations described in Eqs. (gH]) - 
([52|. and the majority of this section will be devoted to 
deriving this system of equations. 

We begin, as in Paper I, with the PN metric at New- 
tonian order, 

dS^ = -(1 - 2UN)df + {1 + 2UN){dR^ + R^d?n) , (24) 

the same as Eq. (g]) of Sec. [Hi though without the gravit- 
omagnetic terms. Here, however, we write the Newtonian 
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potential (expanded to quadrupole order) as 



TT _ 7-r('=0) 



U 



{1=2) 



N 



m=-2 ^ 



(25) 



The first term is the monopole piece (M is the total 
mass of the binary) and the first term in the sum is 
the quadrupole part (and Qm are the quadrupole mo- 
ments of the binary) . These two terms above are identi- 
cal to those of Paper I, but the second term in the sum 
(the polynomial in R with coefficients Fm) is different. 
One can include the terms proportional to Fm, because 
like the Newtonian potential, they are solutions to Pois- 
son's equation. These terms diverge at infinity (which 
restricts their use to the near zone), but they cannot 
be determined from the near-zone dynamics alone, how- 
ever. Burke showed [50| . using the technique of matched 
asymptotic expansions, that the terms with coefficients 
Fm could represent the reaction of the binary in the near- 
zone to radiation losses to infinity. The portion of the 
potential due to the moments f ,„ , therefore, is called the 
Burke- Thorne radiation-reaction potential. 

In the hybrid method, we will find a similar quantity in 
the interior PN spacetime by matching the PN near-zone 
solution to a solution in the Schwarzschild exterior with 
no ingoing waves. Namely, when we assume that there 
are no ingoing waves from past-null infinity in the exte- 
rior BHP spacetime, this determines a radiation-reaction 
potential within the interior PN spacetime. This allows 
us to incorporate the effects of wave propagation in the 
background black-hole spacetime into the dynamics of 
the binary. While the Schwarzschild background does 
not capture every detail of the curvature of a binary at 
small separation, we see that it does capture much of the 
important effects. 

Proceeding with the calculation, we assume we have an 
equal-mass, nonspinning binary in the x-y plane, located 
at 

XA{t) ^-Xsit) = ^ A{t) {cos a{t), sin a{t), 0) , (26) 

where A and B are labels for the two members of the 
binary. Each black hole has mass M/2, and a straight- 
forward calculation shows that 



Q2(t) 

Qo{t) 

Q-2{t) 



37r MAjtf 
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O2W 



(27) 

(28) 
(29) 



where the overline stands for complex conjugate, and 
where the m = ±1 components must be zero for this 
equal-mass binary by symmetry. Throughout this paper, 
we focus just on the m = ±2 multipoles, because as one 
can see from the expressions above, the m = moment 



only evolves due to the radiation-reaction force (for cir- 
cular orbits), and, therefore, is less significant than the 
m = ±2 multipoles, which change on the orbital time 
scale. Moreover, the m = —2 quantity is the complex 
conjugate of the corresponding m = 2 quantity, so when 
we write Q(t) (or any other variable that might be in- 
dexed by m), we refer to the m — 2 variable, and sim- 
ilarly, for Q{t), we mean the m = —2 element. This 
way, the notation can be simplified by dropping the m 
label on multipole coefficients. Thus, we can write the 
quadrupole perturbation as 



U 



N 



u 



N 



Q{t) 

i?3 



where 



Q{t) = 



[3^ MA{tf 



F(t)R^ 



~i2a{t) 



(30) 



(31) 



and F(t), an undetermined function of time, is the 
radiation-reaction potential. 

One can substitute Eq. ((30l) into Eq. (|T9| and use the 
fact that r = i? — A/ to find the Zcrilli function. Cal- 
culating the Zerilli function introduces many factors of 
M/R into the end result, which, because our calculation 
is only accurate to Newtonian order, we will keep only 
the leading-order terms in R. We find that 



2Qit) F{t)R^ 
R^ S 



(32) 



We will also shortly need expressions for the derivative of 
the Zerilli function with respect to the tortoise coordinate 
r*, Eq. ([2|), which we compute here as well. Again, we 
will keep the leading-order expression in R, but we will 
also retain the factor of 



dr 



= l- 



2M\ _ {R- M) 
" {R + M) ' 



(33) 



since although ^(c) rnay be constant on the horizon, 
95'(c)/9r* should vanish there [6l|. The result of this 



calculation is that 



R-M 
R + M 



4QW 

i?3 



F{t)R^ 



(34) 



The Zerilli function satisfies the simple wave equation 
in a potential, Eq. ([21]). As before, the value of the 
Zerilli function at the matching surface, Rs{t) = A{t)/2, 
provides a boundary condition for the Zerilli equation 
on the matching surface, but now there is an additional 
boundary condition on the Zerilli function's derivative 
with respect to the tortoise coordinate. The two bound- 
ary conditions state that 



*(e)(i) 

dr.. 



8Q{t) , Fit)A{t)^ 



A{t)^ 24 
/Ait) - 2M 



\A{t) + 2M^ 

32Q{t) F{t)A{t)^' 



(35) 
(36) 
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By eliminating the unknown function F{t) from the 
above equations, one can impose a mixed (Robin) bound- 
ary condition at the matching surface between the PN 
and BHP spacetimes, 



- 2M 



\A{t) + 2M 



it) 



80Q{t) \ 
A{tr ) 



(37) 

This specifies a boundary condition at a given moment 
in time, but it does not yet describe how to evolve the 
matching surface (through evolving the reduced-mass 
motion of the system) and the value of the Zerilli function 
on this surface. 

One can determine the value of the Zerilli function at 
later times through the boundary condition above, and 
the following additional constraint. By integrating the 
Zerilli equation with respect to u, one finds that 



dv 



+ 



dv 



Vl'-^\r)^^,jiu',v)du' 



(38) 



where we have written r implicitly as a function of u 
and V, and Ugit) denotes the value of u at the matching 
surface for a given time t. Having no ingoing waves forces 
the second term to be zero, so 



dv 



Us(t) 



V,^^-'\r)^^,){u',v)du' . (39) 



Because both 9^'(o)(i)/9w and 9\E'(e)(i)/9r* are con- 
strained at the point of the matching surface, this deter- 
mines the evolution of ^(c)(i) on the matching surface. 
It is easiest to express the Zerilli function on the match- 
ing surface as a function of time via 4'(c)(i, ''*(0)- Then, 
taking the total derivative. 



dt 



^(e)(i) = 



dr^ a* 



dt 



dt 



using the facts that 



9* 



dt 



= 2- 



9* 



(o) 



d-^ 



dv 



dr^ 



and 



It 



dr 
di' 



(40) 



(41) 



(42) 



along with the relationship a(t) 



*(c)(i) 



dv 



2r{t), one can write 
(43) 



lfAitH2M 
2 \A{t) -2M ' ^ ' 



dr^ 



In the above equation, 9^'(e)(t)/9?; is given by the in- 
tegral of the Zerilli function up to that time, Eq. 



and 9^'(o)(<)/5r* is given by the boundary condition, Eq. 
([57| . at that instant. As a result, the only term in Eq. 
that is not yet fixed is the expression for A{t). 
The term A{t) specifies the time evolution of the re- 
duced mass of the binary, which, because it is twice the 
radius of the matching surface between the Schwarzschild 
and PN metrics, could conceivably evolve via either 
the PN equations of motion or those of a particle in 
the Schwarzschild spacetime. We will choose the lat- 
ter, for the same reason as described in Paper I: the 
Schwarzschild Hamiltonian has the advantage that a par- 
ticle falling toward the horizon approaches it exponen- 
tially in time, in the limit that the particle is near the 
horizon. Because we are using this motion to approxi- 
mate the region inside of which PN theory holds, we want 
this space to quickly fall toward the horizon as the theory 
begins to converge slowly. Moreover, the motion should 
move smoothly toward the horizon (so as not to introduce 
high-frequency modes that could escape the black-hole 
effective potential). The PN equations of motion do not 
have these desirable features; we consequently favor the 
point-particle evolution equations in the Schwarzschild 
spacetime. 

We write the evolution equations for the reduced mass 
of the system in their Hamiltonian form. As in Paper I, 
we will describe the dynamics of the reduced mass in PN 
coordinates, because at late times, this causes the point 
particles in the PN metric to approach the horizon in the 
external Schwarzschild spacetime as the reduced mass of 
the system does the same. The equations of motion for 
the reduced mass, /z, are 



A{t) = 

PA{t) = - 



dH 

dpA{t) 
dH 
' dA{t) 



ait) 



dH 



dpa{t) 
Pa{t) = Tait) , 



(44) 



where the Hamiltonian of a point particle in the 
Schwarzschild spacetime is given by 



H{A{t),pA{t),Pc.{t)) 



(45) 



2M\ 



2M\pa(^ PcM 



^ Mt)) 



^i^A{t)^- 



The radiation-reaction force is given by the derivative of 
the radiation-reaction potential with respect to if, and 
it should be evaluated at the location of the matching 
region, 



J'ait) 



1 du'-^=^^'^ 

/i d(p 



'^^Fit)e 



2ia{t)-\ 



where U 



(46) 

represents the quadrupole part of the 
radiation-reaction potential. By solving Eq. ([55]) for F(t) 
in terms of \I'(o)(i) [and because Q{t) is proportional to a 
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real amplitude times e 2iQ(t)^ 



see Eq. (|3T|) ] , one can write 



27r ^(t) 



3[vI/(e)We2-(*)]. (47) 



With the above relationship between the radiation- 
reaction force and the Zerilli function, there is now a 
complete set of evolution equations for the reduced-mass 
motion of the system, the Zerilli function on the matching 
surface, and the Zerilli function in the exterior spacetime. 
This system of equations is given by 



A{t) 

PA{t) 
Pa{t) 



dH 
dpA{t) 
dH 



, ait) 



dH 

dpa{t) 




,2»a(t)l 



(48) 

(49) 
(50) 



{r)'^(c){u ,v)du' 



52^ 



A 



/ 6vl/(e)(t) 80Q(t) \ 



dudv 



(51) 
(52) 



where the Hamiltonian is given by Eq. (j45|) , the potential 
by Eq. ([22]) , and the quadrupole by Eq. (|3T1) . By includ- 
ing a radiation-reaction force, we arrived at a set of evo- 
lution equations that simultaneously evolve the reduced- 
mass motion of the binary and the gravitational waves 
emitted, taking into account the back action of the emit- 
ted radiation on the reduced-mass motion. 



D. Weak-Field Analytical Solution 

First, we will confirm that our procedure recovers the 
correct Burke-Thorne radiation-reaction potential in the 
weak-field limit. If we have an equal-mass binary in a 
circular orbit at a large separation, r* w r i? 3> M , 
then the leading-order behavior of the Zerilli equation, 
Eq. (pij) is just a wave equation in flat space. 



9i2 



^*(e)=0. 



(53) 



If one assumes a product solution (e) = e"^*'(/;(i?), then 
for the radial motion, one must solve the ordinary differ- 
ential equation 



dR2 



' LO -if) 
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(54) 



The solutions for ip/R are spherical Hankel functions 
^/R = h2{ujA/2) = j2{ujA/2) + in2{LuA/2), assuming 



there are no ingoing waves. Here uj corresponds to the 
gravitational-wave frequency. We must match this wave- 
zone solution to the PN near-zone expression for the Zer- 
illi function given by Eq. ([35]); additionally, we must also 
match the derivative of the Hankel function with the ra- 
dial derivative of the PN Zerilli function given in Eq. 
(EH). 

We will write these conditions in the frequency domain, 
where 

B(.)AH,i.Am = (55) 



B{uj)Ah'2{ujA/2) 



A^ 24 ' 

32Q(a;) F{uj)A^ 
l3~ ^ 4 



(56) 



and we must solve for the unknown amplitude B{lij) and 
the radiation-reaction potential F(uj) in terms of the 
quadrupole moment Q{uj) and the spherical Hankel func- 
tion h2{ujA/2). Since the matching takes place at very 
large radii, and, by Kepler's law Acj ~ A~^^^ for circu- 
lar orbits, one can expand the Hankel function in Auj/2. 
This allows one to solve for F as a series in 1/A, whose 
three lowest terms are given by 

If) 8 2 

5l^^'^+25l'''^ + *~^'^ + ^(^"')- 
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The third term is the familiar Burke-Thorne radiation- 
reaction potential (written in the time domain, this is 
proportional to five derivatives of the quadrupole mo- 
ment). The first two terms resemble IPN and 2PN cor- 
rections to the Newtonian potential in the near zone; 
however, these terms represent the effects of time retar- 
dation that are needed to match the near-zone solution 
to an outgoing wave solution in the wave zone. As a re- 
sult, our method recovers, asymptotically, the expected 
result. Consequently, the evolution system, Eqs. - 
(l52t . will also give rise to the correct dynamics in the 
weak-field limit. 



IV. NUMERICAL METHOD AND RESULTS 

We begin this section by describing the numerical 
method that we use to solve the system of evolution equa- 
tions, Eqs. (|48l) ~ (|52p . We then show that the evolution 
equations give rise to reasonable and convergent results. 
With this established, we compare our waveform with 
one from a numerical-relativity simulation, and we close 
this section by interpreting the spacetime of the hybrid 
method. 



Numerical Methods and Consistency Checks of 
the Evolution Equations 



Because the set of evolution equations Eqs. (|48p - ([5^ 
has a somewhat unusual form, we describe our numerical 
method in detail, and we present a few basic checks of 
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discretized grid, it is rare that the Zerilh function on the 
matching surface, denoted by 



FIG. 3: A diagram of how we discretize and evolve the Zerilli 
function. The dots represent the Zerilli function evaluated at 
the grid points (the points of intersection of the dashed lines), 
and the solid black line is the matching surface. For all points 
except those adjacent to the solid line, one can use Eq. (|59p 
directly to numerically evolve the Zerilli function. Near the 
solid line, one can use the same procedure as described in Eq. 
(|59|l . except that one must interpolate the Zerilli function to 
the point "i/s.a to use the same procedure. Further detail is 
given in the text of this section. 



the waveform and its convergence. To find the field out- 
side the matching surface, we use the same method as 
that described in Paper I, a second-order accurate, char- 
acteristic method. If we define the following points on 
the discretized grid (see the portion on the right, away 
from the solid line, in Fig. 

= *{'™(" + Aw, w + Av) , = "^["{{u + Au, v) , 

^E^'f\^;{u,v + Av), ^s^^\:';{u,v), (58) 



then discretizing Eq. ((52)) . one can solve for in terms 
of the other three discretized points and the potential: 



N 



AuAv 

+ *w - *S ^ — Vl^j{rc){^E + ^'iv) 



+0{Au'^Av,AuAv^ 



(59) 



Here Vc is the value of r at the center of the discretized 
grid, {u + Am/2, v + Av/2). 

We must evolve this partial diff'erential equation simul- 
taneously with the five ordinary differential equations 
describing the Zerilli function on the matching surface 
and the surface's position, because all these equations 
are coupled together. We solve the ordinary differential 
equations using a second-order accurate Runge-Kutta 
method. As in Paper I, the Zerilli function along the 
matching surface does not always lie on the uniform grid 
in the u-v plane, and we must be careful when finding 
the Zerilli function at grid points adjacent to the match- 
ing surface. For example, at a given value of u along the 



(60) 



will actually fall along a grid point (see the left side of 
Fig. [3] near the solid line). Similarly, when evolving the 
discretized version of Eqs. - ([5^ . it is again unlikely 
that the Zerilli function along the matching surface at 
the next value of u (advanced by one unit of Au), 

^w,.it + At) = *(c)(u,(i) -f Au,w,(i-H At)), (61) 

will fall at a grid point or even at the same value of v as 
the previous earlier value of the Zerilli function, \l/vi/,s(t). 

To be able to use Eq. (|59p to find the Zerilli function 
at M = Us (t) + Au for the next grid point in v (which we 
denote by \I'jv,s), we must interpolate the Zerilli function 
at fixed u = Us{t) to the same value oi v — Vs{t) as 
+ At). We will label this point by 



^S,s^'i'io}{Usit),V,{t + At)). 



(62) 



As in Paper I, this interpolation does not infiuence the 
convergence of the algorithm when done with cubic inter- 
polating polynomials. With the value of the Zerilli func- 
tion at u = Us(t) and the nearest grid point in v (which 
we will call e,s), one can then find the point ^'^v.s using 
Eq. (1591) . where ^I'w, and are replaced by '^e,s, 
'i'w,s{t + At) and '^s.s, respectively. 

As a final note on the numerical methods, we point out 
that in the evolution equation for the Zerilli function on 
the matching surface, Eq. P5)) . the term d'^ i^c){t)/dv in- 
volves an integral of the Zerilli function times the poten- 
tial, Eq. (131]). Explicitly evaluating this integral adds to 
the computational expense significantly, so we compared 
the value of (c) (t) / dv obtained through performing the 
integral with the value found from evaluating 95'(o) (t) /dv 
numerically using a fourth-order finite-difference approx- 
imation of the derivative, calculated from the Zerilli func- 
tion in the adjacent exterior BHP spacetime. Since the 
two agreed to within the numerical accuracy of our so- 
lution, we used the finite-difference approximation of 
95'(o)(i)/9f in our numerical evolutions. 

We now examine a few consistency checks of the nu- 
merical solutions to the system of evolution equations, 
Eqs. (gS]) - dSll). In Fig. Il we show, in black, the tra- 
jectory of the reduced mass of the binary in the PN co- 
ordinates. On this same figure, we have depicted the 
Schwarzschild black hole by a filled black circle, the light 
ring of this black hole by a red (light) dashed circle, 
and the innermost stable circular orbit (ISCO) by a blue 
(dark) dashed and dotted circle. One can see that the 
radiation-reaction force causes the matching region to 
adiabatically inspiral, until it approaches the ISCO. Once 
at the ISCO, it begins plunging more rapidly toward the 
light ring, and then falls past the light ring and asymp- 
totes to the horizon of the final black hole. 

The initial conditions of this evolution correspond to 
a binary with a PN separation of ^(0) = 14 in a circular 
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FIG. 4: In black, the trajectory of the reduced-mass motion 
of the binary, in the PN coordinate system. The blue (dark) 
dotted and dashed circle shows the Schwarzschild ISCO, and 
the red (light) dashed circle depicts the light ring of the 
Schwarzschild spacetime. The large filled black circle rep- 
resents the horizon. One can see that the binary plunges 
soon after it reaches the ISCO of the exterior Schwarzschild 
spacetime. 



orbit, with no ingoing gravitational waves from past-null 
infinity, and with the radiation-reaction force initially set 
to zero. We do not let the radiation-reaction force en- 
ter into the dynamics (thereby holding the binary at a 
fixed separation) until we have a stable estimate of the 
force. At this point, we include the radiation-reaction 
force (thereby letting the binary begin its inspiral). To 
minimize eccentricity, we introduce a small change in the 
radial momentum ^^(0) that corresponds to the radial 
velocity of a PN binary at that separation. Explicitly, 
we find this value of pA (0) by solving 



i(0) 



dH 



dp Ait) 



t=0 



16 M3 

TA(0)3 



(63) 



(see, e.g. H), while assuming that Pa{0) continues to 
have the value for circular orbits 



Pa(0) = 



MA{0) 



v/A(0)/Af-3 



(64) 



This is necessary to make the orbit as circular as possible 
once the binary begins to inspiral. We do not show the 
initial few orbits before we include the radiation-reaction 
force, and we denote the zero of our time to be the mo- 
ment when we let the radiation-reaction force begin act- 
ing on the binary. 

We also calculate the Zerilli function corresponding to 
these initial conditions, as a function of increasing nu- 



FIG. 5: The norm of the Zerilli function at a given res- 
olution, Av/M, minus the Zerilli function at the highest res- 
olution, {Av)/M = 1/64, which we denote by |^ (c),(Atj)/M — 
'!'((,) i/54|/A'f. We also include a power law proportional to 
[(Ai;)/A/]^ to indicate the second-order convergence of our 
result. 



merical resolution. In Fig. [5l we show that the Zerilli 
function at large constant v, does converge in a way that 
is consistent with the second-order-accurate code we are 
using. We show the norm of the difference between 
the Zerilli function at a given resolution, which we denote 



(c),(A«)/M 



and the highest resolution, {Av)/M = 1/64, 



which we denote by ^'(•, 



we write as hi' 



c),l/64- 



The norm, therefore, 



{c),{Av)/M 



(c), 1/641 



and we normalize 



this by the number of data points in the evolution, and 
the mass. We also include a power law, proportional 
to [{Av)/M]'^, which indicates the roughly second-order 
convergence of the waveform. 

We then plot the real part of the Zerilli function ex- 
tracted at large constant v, for the highest resolution 
{Av)/M = 1/64, in Fig. H The top panel depicts the 
Zerilli function throughout the full evolution. Because 
it is difficult to see the slow increase of the amplitude 
and frequency during early times and the smooth transi- 
tion from inspiral to merger and ringdown at late times, 
we highlight the early stages of the inspiral in the lower- 
left panel, and we depict the merger and ringdown in 
the lower-right panel. Because \/6^(c) = — ihy,), 

for the I = 2 modes at large r [see Eq. (jlOO|) ]. the Zerilli 
function is essentially identical to the gravitational wave- 
form. From this one can see the hybrid method produces 
a smooth inspiral-merger-ringdown waveform. Because 
the hybrid waveform has the correct qualitative features 
of a full inspiral-merger-ringdown waveform, it is natu- 
ral to ask how well it could match a numerical-relativity 
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FIG. 6: The top panel shows the real part of the Zerilli func- 
tion throughout the entire evolution, extracted at large con- 
stant V. The bottom- left panel displays the early part of the 
same Zerilli function, and the bottom-right zooms in to the 
merger and ringdown portions of the function. Because only 
a factor of \/6 differentiates the Zerilli function from r times 
the waveform, this can be thought of as the waveform as well. 



waveform. We, therefore, turn to this question in the 
next section. 



B. Comparison with Numerical Relativity 

In this section, we will first discuss how well the wave- 
form compares with a similar waveform from numerical- 
relativity simulations. The first part of the section is 
devoted to showing how we can make small modifica- 
tions to the hybrid procedure to make the phase agree 
well with that of a numerical-relativity waveform dur- 
ing inspiral (though the comparison of the amplitudes is 
less favorable). The second part of this section describes 
why the hybrid method, in its current implementation, 
does not agree well with numerical-relativity simulations 
during the merger and ringdown phases. The reason for 
the discrepancy during the late stages of the waveform is 
well understood (the background spacetime of the hybrid 
method is Schwarzschild, whereas the final spacetime of 
the numerical simulation is Kerr) and could be improved 
by modifications to the hybrid method. 



1. Agreement of the Waveforms during Inspiral 

We will briefly describe a small change to the hybrid 
method that leads to a waveform whose phase agrees 



well with a numerical-relativity waveform during the in- 
spiral part. We will continue to find the Zerilli func- 
tion through the procedure describe in Sec. IIIICI using 
the leading-order expression for the Newtonian potential 
(and thus also the leading-order radiation reaction). We 
note, however, that when we took the derivative of the 
Zerilli function on the matching surface with respect to 
Eq. dMl), we kept the factor of (1 - 2M/r). This is 
reasonable, physically, because, although the Zerilli func- 
tion itself may approach a constant on the horizon, its 
derivative with respect to should vanish. Conversely, if 
the derivative of the Zerilli function did not vanish, then 
that could correspond with a perturbation that diverges 
on the horizon. Nevertheless, because the boundary con- 
dition only takes into account the leading Newtonian ex- 
pressions, the overall factor of (1 — 2M/r) is a higher 
PN correction, from the point of view of the interior PN 
spacetime. We, therefore, are justified in dropping this 
term in our leading Newtonian treatment, and we find 
the agreement between numerical relativity and the hy- 
brid method is helped by this. It is likely that further 
adjustments will lead to even better results, though a sys- 
tematic study of this is beyond the scope of this initial 
exposition. 

The modification above results in only a small change 
to Eq. (IMl), 



d^^,){t) _ 32Q(t) F{t)A{ty 



A{tr 



(65) 



and it also alters the boundary condition, Eq. (|37)) of Sec. 



6 



A{t) 



*(c)(^) 



8og(t) 

A{tr 



(66) 



With the exception of these two equations and the fact 
that we begin the evolution from a larger initial radius, 
A{0) = 15.4, we evolve the new system of equations in 
exactly the same way as that described in detail in Sec. 
CVAl 

For our comparison with a numerical-relativity wave- 
form, we use the I = 2, m = 2, mode of the waveform 
from an equal-mass, nonspinning, black-hole binary de- 
scribed in the paper by Buonanno et al. fl8|. In this 
simulation, the black holes undergo 16 orbits before they 
merge, and the final black hole rings down. We plot 
the numerical-relativity waveform in black in Fig. [71 and 
we show the equivalent waveform from our approximate 
method in red (gray). Recall that the I = 2 modes of the 
Zerilli function are related to the waveform by 



(67) 



[see Eq. dlOOp ]. Although the amplitudes of the wave- 
forms do not agree exactly, the fact that the phases match 
so well throughout the entire inspiral is noteworthy. The 
approximate waveform completes one more orbit than 
the numerical-relativity one, and the ringdown portions 
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FIG. 7: In black is the real part of the I — 2, m = 2 mode 
of a numerical-relativity waveform, whereas in red (gray) is 
the equivalent quantity from the approximate method of this 
paper. The agreement of the waveforms' phases is quite good 
throughout the entire inspiral, although the amplitudes dif- 
fer. The approximate and numerical-relativity waveforms dif- 
fer during ringdown, because the approximate method uses a 
black-hole with zero spin, whereas the final black hole in the 
numerical-relativity simulation has considerable spin. 



differ as well. This is not too surprising, however, since 
the final black hole in the numerical-relativity simulation 
is a Kerr black hole with dimensionless spin x ~ 0.7 (see, 
e.g., Scheel et al. [51]), whereas our ringdown takes place 
around a Schwarzschild (nonspinning) black hole. 



2. Differences in the Instantaneous Frequency during 
Merger and Ringdown 

The discrepancy between the two waveforms at late 
times in Fig. [7] is most evident in the instantaneous fre- 
quency, often defined as 



parable to this constant, there are spurious oscillations in 
the frequency as it becomes dominated by this constant 
offset. The hybrid waveform needed no modification. 

Solid curves depict the instantaneous frequency of the 
numerical-relativity waveform in Fig. [H] the real (oscil- 
latory) part is the black curve and the imaginary (de- 
caying) part is the red (gray) curve. Similarly, the 
black dashed curve is the real part of the instantaneous 
frequency of the hybrid method, and the red (gray) 
dashed curve is its imaginary part. The hybrid and the 
numerical-relativity frequencies are in very good agree- 
ment for the inspiral up until the late stages highlighted 
here. The numerical-relativity waveform quickly tran- 
sitions after the plunge and merger to the least-damped 
I = 2^ m — 2 quasinormal-mode frequency and decay rate 
for a Kerr black hole of final dimensionless spin equal to 
roughly x ~ 0-7 (see, e.g., 64]). The frequency of the 
hybrid waveform, however, undergoes a similar qualita- 
tive transition, but it approaches the least-damped ? = 2, 
m = 2 ringdown frequency of a non-spinning black hole 
(the background of the hybrid method). The hybrid 
method, however, oscillates around this value with a fre- 
quency that is proportional to twice the frequency of this 
least-damped, I = 2, m = 2 quasinormal mode. 

The origin of this oscillation is simple and, in fact, was 
explained by Damour and Nagar [38| . For each I and 
TO, there are quasinormal modes with both positive and 
negative real parts, which both have a negative decay 
rate. For a Schwarzschild black hole, the decay rates 
are the same and the real frequencies are identical, but 
have the opposite sign. For a Kerr black hole, however, 
the positive-frequency modes have a lower decay rate 
than the negative-frequency modes (and the positive fre- 
quency is larger in absolute value than the negative fre- 
quency is). While a counter-clockwise orbit will tend to 
excite predominantly the mode with a positive real part, 
it can also generate the negative real-frequency mode as 
well. In the hybrid waveform, because the background 
is Schwarzschild, the positive- and negative-frequency 
modes decay at the same rate, and they can interfere 
to make the oscillations at twice the positive real fre- 
quency. In the numerical-relativity waveform, however, 
the difference of the frequencies and decay rates prevents 
this from happening. 



Mu 



(68) 



where "^(c) is the Zerilli function measured at large r. 
We calculate this frequency for both the hybrid and the 
numerical-relativity waveforms, and we show the real and 
the imaginary parts (the oscillatory and damping por- 
tions, respectively) in Fig. [S] The numerical-relativity 
waveform was offset from zero at late times by a small 
constant of order 10~^. We subtracted this constant from 
the waveform to find the instantaneous frequency; other- 
wise, when the amplitude of the waveform becomes com- 



C. Interpreting the Hybrid Waveform and 
Spacetime 

Since the phase during inspiral agrees so well, and be- 
cause the transition from inspiral to merger and ring- 
down is qualitatively similar, this leads one to wonder 
to what extent the hybrid approach may also be a use- 
ful tool for generating gravitational-wave templates for 
gravitational-wave searches. To capture the correct ring- 
down behavior, the hybrid method would need to be 
extended to a Kerr background; however, it is likely 
that calibrated approaches using the effective-one-body 
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FIG. 8: The solid curves are the instantaneous frequency 
[see Eq. H68|) ] of the numerical-relativity waveform; the black 
curve is the real, oscillatory part and the red (gray) curve 
is the imaginary, decaying part. The black, dashed curve 
and the red (gray) dashed curve are the real and imaginary 
parts, respectively, of the frequency for the hybrid method. 
The frequencies agree quite well during the inspiral, but at 
late times they begin to differ. The qualitative transition 
from inspiral to merger and ringdown is similar, but the final 
quasinormal-mode frequencies that the waveforms approach 
differ, because the numerical-relativity simulation results in 
a Kerr black hole of dimensionless spin x ~ 0.7, whereas the 
hybrid waveform is generated on a Schwarzschild background. 
The oscillations in the hybrid waveform arise from the inter- 
ference of positive- and negative-frequency modes that can 
arise in a Schwarzschild background, as explained in the text 
of this section. 

method (see, e.g., [ll]) or phenomenological frequency- 
based templates (see, e.g., [631) "^ill be more efficient for 
these purposes. The hybrid approach, as described here, 
will likely be more helpful as a model of how the near- 
zone motion of the binary connects to different portions 
of the gravitational waveform. 

As an example of this, we show the real part of the 
gravitational waveform at large v, the black solid curve, 
and the corresponding value of the Zerilli function on the 
matching surface, the red (gray) dashed curve in Fig. [HI 
Interestingly, the Zerilli function on the matching surface 
and that extracted at large constant v are roughly out-of- 
phase with one another during the inspiral; namely, along 
a ray of constant u, the Zerilli function undergoes nearly 
one half cycle as it propagates out to infinity. This fea- 
ture is also visible in Fig. [TOl but it is harder to discern 
there. This behavior holds through inspiral up to the 
beginning of the merger. During the merger, however, 
the two transition away from the out-of-phase relation- 



FIG. 9: The real part of the Zerilli function on the matching 
surface, the red (gray) dashed curve, and the real part of 
the gravitational waveform, proportional to the real part of 
the Zerilli function at large v, (the black solid curve). The 
two functions are nearly out-of-phase for the inspiral, and 
the wave propagates more or less directly out. During the 
merger, they begin to lose this phase relationship, and during 
ringdown the Zerilli function on the matching surface becomes 
constant. This implies that the ringdown waveform is due just 
to the waves scattered from the potential, as also illustrated 
in Fig. [To] 



ship, before the Zerilli function on the matching surface 
becomes a constant during the ringdown (when the re- 
duced mass of the binary falls toward the horizon along 
a line of constant v). 

This change in phasing between the Zerilli function on 
the matching surface and that at large v (along a line 
of constant u) allows one to give an interpretation to 
the different parts of the waveform. The inspiral occurs 
when the waveform propagates out directly, but nearly 
out-of-phase with the matching surface. The merger is 
the smooth, but brief, transition during which the phase 
relationship between the matching surface and the wave- 
form evolves, and the ringdown is the last set of waves 
that are disconnected from the behavior on the surface 
(they are the scattered waves from the potential barrier) . 

We also show in Fig. [TU] a contour-density plot of the 
real part of the Zerilli function in the u-v plane during 
the last few orbits of inspiral, the merger, and the ring- 
down (for the evolution discussed in this section). This 
is a spacetime diagram, where time runs up, and the ra- 
dial coordinate, increases to the right. The matching 
surface is the dark timelike curve running up that turns 
to a line of constant v at the end. The region to the 
left of the surface, the solid green (gray) is the interior 
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FIG. 10: A contour-density plot of the real part of the Zerilli function for the evolution discussed in this section. We only show 
the last few orbits of the inspiral, followed by the merger and ringdown. In this spacetime diagram, time runs up, r, increases 
to the right, and the coordinates u and v run at 45 degree angles to the two. The line that starts at nearly constant t and 
evolves to a line of constant v is the matching surface, and to the left of this line, the solid green (gray) region is the interior PN 
region (where we do not show any metric perturbations). The exterior is the BHP region, where we show the Zerilli function. 
During inspiral, the Zerilli function propagates out almost directly, and it oscillates between positive, yellow (light gray) colors, 
and negative, light blue (darker gray). Black dashed contour curves are used to highlight this oscillation. As the reduced mass 
of the binary plunges into the potential during merger, the amplitude and frequency of the radiation increases, but it promptly 
rings down to emit little radiation, in the upper green (gray) diamond of the diagram. There are black dashed contours here as 
well to indicate that there is still oscillation, even though it is exponentially decaying (and hard to see through the color scale). 



PN region, but we do not show the metric perturbation 
in this region. On its right is the BHP region, where 
we show the Zerilli function colored so that blue colors 
(dark gray) are negative and red colors (light gray) are 
positive. Away from the matching surface, the Zerilli 
function oscillates between yellow (light gray) and light 
blue (darker gray) for several orbits before inspiral. Each 
oscillation is bounded between a black, dashed contour 
curve. As the reduced mass of the binary plunges toward 
the horizon, the outgoing waves increase in frequency and 
amplitude, which is how we describe the transition from 
the inspiral to the merger phase. The merger phase is 
short, and the black hole rings down (leading to very lit- 
tle gravitational- wave emission in the top corner of the 
diagram) . As the reduced mass of the system approaches 
the horizon, there is a small wavepacket of ingoing radi- 
ation that accompanies it. 

We close this section with one last observation. If we 



were to plot the equivalent quantities to those in Figs. [9] 
and [TU] for the evolution in Sec. Illli then one would see 
that the Zerilli function on the matching surface increases 
during ringdown instead of approaching a constant. This 
does not have any effect on the waveform, because it is 
a low frequency change that occurs within the potential 
barrier, and is hidden from the region of space outside the 
potential. In some sense, it is a strong confirmation of 
Price's idea that the details of the collapse will be hidden 
within the potential barrier. At the same time, however, 
this behavior arises from the fact that the derivative of 
the Zerilli function with respect to vanishes on the 
matching surface. When this condition was neglected in 
this section, it led to a more regular behavior there. This 
suggests that it may be worth while to do a more careful 
analysis of how the Zerilli function and its derivatives 
near the horizon should scale in the presence of radiation 
reaction. 
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V. SPINNING BLACK HOLES, SPIN 
PRECESSION, AND THE SUPERKICK MERGER 

In this section, we will incorporate the effects of black- 
hole spins into our method, with the aim of understand- 
ing the large kick that arises from the merger of equal- 
mass black holes with spins antialigned and in the orbital 
plane (the superkick configuration). To do this, we will 
first discuss adding odd-parity metric perturbations to 
the results in the previous section. We will then indicate 
why spin precession is important in producing large kicks 
and discuss two ways of implementing spin precession: 
the PN equations of precession and geodetic precession 
in the Schwarzschild spacetime. In our method, we will 
use the geodetic-precession approach, and we will present 
numerical results for the kick that uses this equation of 
spin precession. 



Odd-Parity Metric Perturbations of Spinning 
Black Holes 



To incorporate the effects of spin into our model, we 
will add the lowest-order metric perturbation arising 
from using spinning bodies in the PN metric, as we did 
in Paper I. This comes from the metric coefficients 



hoi — 



(69) 



Here we use the notation of Paper I, where we label the 
two bodies by A and B. The new variables S\ represent 
the spin angular momentum of the body, Ra is the dis- 
tance from body A and n\ is a unit vector pointing from 
body A. The variables for body B are labeled equiva- 
lently. Since we will focus on the extreme kick configura- 
tion, we will assume the black holes lie in the x-y plane, 
at positions Xyi(f) and ^B{t) [identical to Eq. (|26|) of 
Sec, nil C] . and that the spins are given by 

SA{t) = -SBit) = 5(cos/3(0,sin/3(t),0) , (70) 

where S — x(M/2)^ is the magnitude of the spin, and x 
is the dimensionless spin, ranging from zero to one. 

Under these assumptions, one can show that the Carte- 
sian components of the metric coefficients above are 

hox = -^^^^ sm29 sin l3{t) cos{a{t) - (p) , (71) 
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hoz = 
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sin 29 cos (3{t) cos{a{t) - (p) , (72) 
sin[a{t) - I3{t)] (73) 
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spherical-polar coordinates to find that 
2SA{t) 



h-OR — 



hoe — — 



i?3 

2SA{t) 

i?2 



sin[a(i) — f3{t)] cos 9 , 
sin[a{t) ~ I3{t)] cos 9 
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(74) 
(75) 

), 

(76) 



As written above, the metric perturbations do not take 
the form of an odd-parity vector harmonic, because there 
is a dipole-like piece in two of the components. This can 
be eliminated by making a gauge transformation, 

Co = -^^cos0sin[a(t)-/3(t)]. (77) 

A small gauge transformation produces a change in the 
metric via 



^^^u J 



(78) 



which in this case sets /lofl — 0. The remaining terms 
in the metric can then be expressed in terms of the odd- 
parity, vector spherical harmonics. 



-2, ±2 
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(79) 



X^'O = ix','",X^/) = -^\l^sm'9cos9{0^)m 



2 V TT 



A short calculation shows that 



{hoe, h 




8SA{t) /!Lcos[a(t)-/3(i)]X2'0 



One can then convert the Cartesian components into 



As with the even-parity, mass-quadrupole perturba- 
tions discussed in the previous section, we will only be 
interested in evolving the m = 2 perturbation (though 
in this case it is an odd-parity, current-quadrupole mo- 
ment). The reason for this is subtle, and will be clar- 
ified in the next section. Nevertheless, we will mention 
here that during the merger and ringdown (when the kick 
is generated), the spins precess at the orbital frequency 
[namely a{t) = ${t)]. As a result, the m = part of 
the perturbations which depend on a{t) — f3{t) become 
constant, and the only changes in the perturbations come 
from changes in A{t). We mentioned in Sec. lIII Cl that we 
would also neglect the m — part of the even-parity per- 
turbations, because it also evolved from time variations 
in A{t)^ which occur on the time scale of the radiation- 
reaction force (2.5 PN orders below the leading-order or- 
bital motion). Consequently, because we are interested 
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in the behavior of the binary during merger and ring- 
down, we can neglect the m — parts of the odd-parity 
metric perturbations for this same reason. In addition, 
because we are treating just the m — ±2 perturbations 
(and the m = —2 term is the complex conjugate of the 
TO = 2 moment), we will again drop the label m on the 
perturbations. 

Thus, the relevant piece of the gravitomagnetic poten- 
tial for our calculation will be 



w 



(o) 



SA{t) 



^g-»[a(t)+/5(t)] 



(82) 



and one can then use Eq. (|20p and the fact that R — r~M 
to find that the Regge- Wheeler function is (at leading 
order in r), 
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-ila(t)+l3{t)] 
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This means that on the matching surface. 
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(84) 



We can then evolve the Regge- Wheeler equation, Eq. 
(PT|) . (with the odd-parity I = 2 potential) using Eq. 
([84)) as the boundary condition along the matching sur- 
face. We will not take any radiation-reaction effects from 
the current-quadrupole perturbations into account (since 
they are 1.5 PN orders below the leading-order Newto- 
nian radiation reaction of Sec. IIII CP : as a result, we will 
evolve the Regge- Wheeler function using the matching 
surface generated by the even-parity, mass-quadrupole 
perturbations alone. 



B. Spin Precession 

Before we discuss the evolution of the Regge- Wheeler 
and Zerilli functions, we will mention an effect that is 
important for our recovering the correct qualitative be- 
havior of the kick in superkick simulations. This effect 
was observed by Schnittman et al. in [s^l and clarified to 
us by Thorne [55|. In Schnittman et al.'s discussion of 
the superkick configuration, the authors observe that the 
spins precess in the orbital plane very rapidly during the 
merger, approaching the orbital frequency just before the 
ringdown. We will give a heuristic argument of why this 
effect should occur before we explore two models that 
produce spin precession (one based on the PN equations 
of motion and the other based on geodetic precession in 
the Schwarzschild spacetime). We will ultimately favor 
the latter. 



1. Motivation for Spin Precession 

One can see the need for spin precession from the fol- 
lowing simple argument. Just as the even-parity pertur- 
bations gave rise to a waveform that increased from twice 



the orbital frequency to the quasinormal-mode frequency 
during the merger phase (see Fig. [5]), so too must the odd- 
parity perturbations of the previous section give rise to a 
part of the waveform that transitions from the orbital fre- 
quency to the same quasinormal-mode frequency as the 
even-parity perturbations. The quasinormal-mode fre- 
quencies are the same, because both the Regge- Wheeler 
and Zerilli functions are generated by Z = 2, m = ±2 
perturbations. Because the Zerilli function is generated 
by a boundary condition proportional to e"*^"'^*^ and the 
Regge- Wheeler function produced by a boundary condi- 
tion that changes as e"*!"^*-*"*"^^'^!, for the two perturba- 
tions to evolve in the same way, both a(t), the orbital 
evolution, and /3(i), the spin precession, should evolve in 
identical ways at the end of merger. Stated more phys- 
ically, at the end of merger, the spins should precess at 
the orbital frequency. 

This rapid precession of the spins was observed by 
Briigmann et al. (56j in their study of black-hole super- 
kicks. Using a combination of PN spin precession and 
numerical-relativity data, they were able to match the 
precession of the spin in their numerical simulations. We 
will explain in the next section why this worked so well 
for their simulation, but why it will not work as well in 
the hybrid method. 



2. Post-Newtonian Spin Precession 

Briigmann et al. begin from the well-known spin pre- 
cession for a binary (see, e.g., [66j). 

^^(^^ = + ^) f^^^^^ ^^^^^^ ^^^^ 

where we just write the leading-order effect from the 
Newtonian angular momentum, 

L^it) = M[Xa(0 - XB(i)] X [±A{t) - Xb(0]} (86) 

The vector n is a unit vector from the center of mass. 
There is an equivalent equation for the precession of 
SB(t), identical to the equation above, under the in- 
terchange of A and B. Given the form of the equation 
above, the magnitude of the spin does not change, and 
the spin precesses about the Newtonian angular momen- 
tum LAr(i). Moreover, Briigmann et al. found that for 
the superkick configuration, where the spins lie in the 
plane, precession of the spins does not produce a large 
component out of the plane (the z component in this 
case). 

For simplicity, therefore, we will just consider the com- 
ponents of the spin in the orbital plane, which, at leading- 
order, will precess as a result of coupling to the Newto- 
nian orbital angular momentum. The Newtonian angular 
momentum is 



LAr(i) = ^iA{t)^a{t)z, 



(87) 
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where a{t) is the orbital frequency. With the assumption 
that S\ = Sg = 0, the spins precess via the equation 



7Ma{t) 
8A{t) 



z X SA{t)] 



where we have also used the fact that this is an equal- 
mass binary, {Ma = Mb = M/2 and n = M/4). Taking 
the time derivative of Eq. (j70|) . we obtain the expression 
for the left-hand side of the equation above, 



SA{t) ^ X SA{t)] 



(89) 



Relating the two expressions, we arrive at the equation 
of spin precession. 



7M 

8A{t) 



a(t) . 



(90) 



For the hybrid method, this expression will not lead to 
the spin-precession frequency approaching the orbital fre- 
quency, since A{t) > 2M for the entire evolution (and 
hence, the spin-precession frequency will not even be 
half the orbital frequency at its maximum). In the 
next section, we will put forward an equation of spin 
precession based on geodetic precession in the exter- 
nal Schwarzschild spacetime, which will have the desired 
spin-precession behavior. 

Before turning to the next section, we address the ques- 
tion of why PN spin precession worked so successfully for 
Briigmann et al. Their initial data begins in a gauge that 
is identical to the 2PN ADMTT gauge, and they assume 
that it continues to stay in that gauge throughout their 
evolution. As a result, they use the puncture trajectories 
as the positions of the black holes, and the 2PN ADMTT 
gauge expressions to relate the momenta of the black 
holes to their velocities. Although the PN equations of 
spin precession are written in harmonic gauge, they use 
the puncture results to calculate these expressions. This 
is reasonable, because the harmonic and ADMTT gauge 
positions do not differ much until separations of roughly 
A{t) « 2M. Their puncture separations do reach small 
values of A{t) < M prior to merger, and they continue 
to use the harmonic-gauge spin-precession formula in this 
regime (even as the PN approximation starts becoming 
less accurate). This works remarkably well, nevertheless, 
and, as one can see from Eq. (|90p . when A{t) ss 7M/8, 
the spins will precess at the orbital frequency. Thus, the 
work of Briigmann et al. helps to confirm that the locking 
of the orbital and spin-precession frequencies is impor- 
tant in the superkick merger, but to replicate this effect 
in the hybrid method will require a different approach, 
described below. 



3. Geodetic Precession in a Schwarzschild Spacetime 

Our approach to spin precession relies on geodetic pre- 
cession in the Schwarzschild spacetime, which we review 



below. The problem of geodetic precession of a spin on 
a circular orbit in the Schwarzschild spacetime is well 
understood; its derivation appears in the introductory 
text by Hartle [i^l, for example. We will reproduce 
some of the important elements of the derivation here, 
using our notation, however. One typically starts with 
the spin 4-vector S'^ (whose spatial components lie in 
the orbital plane) that travels along a circular geodesic 
parametrized by a 4- velocity u^^. As usual u^u^ = — 1, 
and one also imposes the spin-supplementary condition, 
= 0. The components of these two vectors are 
S = (5^S'^0,S"^), and u = u*(l, 0, 0, d(t)). Because of 
the spin-supplementary condition and the normalization 
of the four velocity, the components S** and u* are not in- 
dependent variables. Thus, when one writes the equation 
of geodetic precession of the spin [Eq. (14.6) of Hartle], 



dT 



(91) 



for circular equatorial orbits, it reduces to two coupled 
equations for the independent variables S'^{t) and S"^(t) 
[Eqs. (14.3a) and (14.3b) of Hartle], 



S"'(i) 



[rs{t) ~?,M]a{t)S'^{t) = 0, (92) 

S^[t) + ^S'-{t) = 0. (93) 
rs{t) 



The dot still refers to derivatives with respect to coordi- 
nate time t (not proper time r). If we assume that a{t) 
does not change much over an orbit (which is true dur- 
ing most of the evolution of the binary, as it changes only 
due to the radiation- reaction force), and we continue to 
denote the angle of the spin in the orbital plane by /3(i), 
then one can write the solution to these equations [Eqs. 
(14.16a) and (14.16b) of Hartle] as. 



/ 2 M 

S'-it) = sJl--^cos[a{t)~m], (94) 



S 1^ 2M d(t) 



rs{t)\] rs(t) a{t) - P{t) 
X sm[a{t) - (i{t)] , (95) 



where the spin is normalized S^S^i, — S'^, and where [Eq. 
(14.15) of Hartle] 



d(t)-/?(i) - Jl- 



3M 



(96) 



Because we only describe the spins with leading-order 
physics, we will only keep the leading-order behavior of 
the spins. Thus, we will describe the spatial components 
of the spins by 



S''{t) ^ Scos[a{t)-P{t)], S'^it) 



S 



rsit) 



sm[a{t)-l3{t)] , 
(97) 
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and we will expand the equation for the evolution of /3{t) 
in a Taylor series up to linear order in M/rg{t), 



1 



3A/ 



We ultimately arrive at the expression that we will use 
to describe spin precession in our formalism, 



Ml 

m 



a{t), 



(99) 



since at leading order A{t) = a{t). 

Although Eq. ([M)) looks quite similar to the leading- 
order PN spin precession, Eq. (j90p . the former equation 
produces a much stronger spin precession than the lat- 
ter does. Not even the next-order PN spin-precession 
terms will produce such strong precession (see, e.g., [68j). 
The equation of spin precession based on geodetic mo- 
tion takes on more of the strong-gravity character of the 
Schwarzschild spacetime. It states that when a spinning 
particle orbits near the light ring, its spin will lock to the 
its orbital motion. An effect quite similar to this happens 
during the merger phase in the superkick simulation, as 
was shown in the work of Briigmann et al., and which 
we discussed in the previous section. In the next section, 
we will show how this contributes to the large kick of the 
superkick simulations. 



Numerical Results and Kick 
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FIG. 11: The top panel shows the real part of the Regge- 
Wheeler function throughout the entire evolution. The 
bottom-left panel just focuses on the early times of inspi- 
ral, where the Regge- Wheeler function slowly increases in fre- 
quency and amplitude because of the binary's inspiral and the 
slow spin precession. In the bottom-right panel, one sees that 
as the spins begin to precess near the orbital frequency, the 
Regge- Wheeler function dramatically increases in amplitude 
and frequency. 



In the first part of this section, we describe how we nu- 
merically solve the Regge- Wheeler equation (we continue 
to solve the Zerilli equation in the same way as described 
in Sec. IIV[) . and we show a representative waveform ob- 
tained from the Regge- Wheeler function. We next de- 
scribe how we calculate the linear-momentum flux and 
the kick from the waveforms. Finally, we close this sec- 
tion by studying the dependence of the kick on the initial 
angle between the spins and the linear momentum of the 
PN point particles. We recover results that are quali- 
tatively similar to those seen in full numerical-relativity 
simulations. 



1. Numerical Methods and Waveforms 

To calculate the Regge- Wheeler function, and thus the 
radiated energy-momentum in the gravitational waves, 
we first make the following observation. Because the odd- 
parity perturbation of the spins of the black holes is a 1.5 
PN effect, the corresponding radiation-reaction force will 
also enter at 1.5 PN beyond the leading-order radiation- 
reaction force discussed in Sec. IIII CI Consequently, we 
do not take it into account in the leading-order treat- 
ment of the radiation-reaction force. Moreover, we note 
that the spin-precession angle, /3(i), does not enter into 
the evolution equations for the reduced mass or for the 



Zerilli function. As a result, the evolution of P{t) and 
^'(o) can be performed after the evolution of the binary 
without spin. In fact, the evolution of (o) is carried out 
in the same manner as that described in Paper I, because 
the matching surface is driven by radiation-reaction from 
the even-parity Zerilli function alone. Were we to include 
the radiation reaction arising from the spins, however, we 
would need to evolve the equations for /3(t) and '^{o) si- 
multaneously, and in a manner identical to that described 

in Sec. nira 

Our initial conditions are identical to those described 
in Sec. lIV Al but we will set the dimensionless spin x = 1, 
and let /3(0) vary over several values from to 27r, to 
study the influence of the initial angle on the kick. We 
first show the real part of the Regge- Wheeler function 
extracted at large constant w, in Fig. [TTJ The top panel 
is the full Regge- Wheeler function, whereas the bottom- 
left panel features the early part from the inspiral (so 
that one can see the gradual increase in the amplitude 
and frequency that comes from the combined effects of 
the binary inspiral, and the increased rate of spin pre- 
cession). In the bottom-right panel, we show the merger 
and ringdown phase, which is obscured in the top panel. 
As the spins start precessing near the orbital frequency 
during merger, one can see the rapid growth of the Regge- 
Wheeler function. 

To see how this spin precession leads to a large kick. 
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FIG. 12: The top panel shows the real part of the Zerilli 
function, the red (gray) curve, and the imaginary part of the 
Regge-Wheeler function, the black curve, throughout the in- 
spiral, merger, and ringdown. The bottom-left panel shows 
only the inspiral, where the Regge-Wheeler function is much 
smaller than the Zerilli function, and oscillates at approxi- 
mately half the frequency. In the bottom-right panel, during 
the merger and ringdown, as the spins precess near the orbital 
frequency, the Regge-Wheeler function increases in amplitude 
and frequency, and becomes in phase with the Zerilli function. 
This leads to a large kick. 



we plot both the even- and the odd-parity metric pertur- 
bations extracted at large constant v in Fig. [T^l We show 
the real part of the Zerilli function, ^(c), in red (gray) 
and the imaginary part of the Regge-Wheeler function, 
\l/(o), in black, for the angle /3(0) that gives the maximum 
kick. As we show below, in Eq. (|103|) . it is the relative 
phase of the product of these components that is impor- 
tant in producing the kick. During the early part of the 
evolution, the Regge-Wheeler function is quite small and 
oscillates with roughly half the period of the Zerilli func- 
tion. This is difficult to see in the upper panel of the full 
waveforms in Fig. 1121 but is more evident in the lower-left 
panel, showing just the early parts of the evolution. In 
the last orbit before the merger and ringdown (shown in 
the lower-right panel), the spins start precessing rapidly, 
and, in the case that produces the maximum kick, the 
real part of the even-parity perturbation function, and 
the imaginary part of the odd-parity function oscillate 
in phase during the merger and ringdown. (For the case 
with zero kick, the two functions are now out-of-phase 
by 90 degrees.) 



2. Calculation of the Kick 

We now discuss, more concretely, how we calculate 
the kick emitted in gravitational waves. At radii much 
larger than the reduced gravitational wavelength, r 3> 
AGw/(27r), one can relate the gravitational-wave polar- 
izations h+ and to the Regge-Wheeler and Zerilli 
functions via the expression. 



2r 



i{i + 2y. 

(1-2)1 



l,m 



-2Ylm , (100) 



where -2i^m is a spin-weighted spherical harmonic. The 
energy radiated in gravitational waves is typically ex- 
pressed as 



lim 

r->-oo 167r 



ihyl^dn. 



(101) 



where is a radial unit vector and dVl is the area element 
on a 2-sphere. A somewhat lengthy calculation can then 
show that the momentum flux in the z direction is given 
by 



PAt) = 



1 



IGtt ^ 2(1 



E 



a + 2)! 



2)! 



i^/,m^/ + l,m 
(o) (o) 



where ci^m — 2m/[Z(/+l)], and di^m is a constant that also 
depends upon I and m. The equations above appear in 
several sources; these agree with those of Ruiz et al. [58 1 
[see their Eqs. (84), (11), (94), and (43), respectively]. 

In our case, however, we just treat the I = 2 and 
TO = ±2 modes of the Regge-Wheeler and Zerilli func- 
tions, and the momentum flux coming from these modes 
greatly simplifles. Because the to = ±2 modes are com- 
plex conjugates of one another, we flnd that the momen- 
tum flux is 



1 



P,(t) = -5[^'(e)*(o) 



(103) 



When we discuss the kick velocity as a function of time, 
we mean that we take minus the time integral of the 
momentum flux, normalized by the total mass, i.e. 



kick 



1 



Pz{t')dt' . 



(104) 



We continue to normalize by the total mass M, be- 
cause numerical-relativity simulations have shown that it 
changes only by roughly 4% durin g a black-hole-binary 
merger (see, e.g., Campanelli et al. [52l|): as a result, nor- 
malizing by the total mass M will not be a large source 
of error. 
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FIG. 13: The momentum flux, Pz{t), as function of time, for 
several values of /? — /3o, where /3o = 2157r/192 is the value 
that gives zero kick. We also include a straight green (light 
gray) dashed line at zero flux to indicate how the momentum 
flux varies around this point. 



3. Numerical Results of the Momentum Flux and Kick 

We now show the results of our numerical solutions for 
the superkick configuration. We first show in Fig. [T^ the 
momentum flux for several different initial angles of the 
spins, (3. In the plots, we subtract the value that gives 
zero kick, which we denote by /3o = 2157r/192. While 
the shape of the pulse of momentum flux has a similar 
shape to that seen in numerical-relativity simulations by 
Briigmann et al. [sgI ]. the absolute magnitude is some- 
what larger. 

The increased overall magnitude of the kick becomes 
more apparent when we plot i'kick(0 Fi§- 1141 where 
■fkickCi) is deflned by Eq. (jl04l) . As one can see, the 
largest value of the kick is near 0.08 in dimensionless 
units, which is roughly a factor of 6 times larger than 
the estimated maximum from numerical-relativity sim- 
ulations at lower dimensionless spin parameters. This 
is largely because the even-parity Zerilli function (pro- 
portional to the waveform) is also signiflcantly larger in 
amplitude than that of numerical-relativity simulations. 

Nevertheless, we then show, in this model, that the 
kick depends sinusoidally upon the initial orientation 
of the spins, as seen in numerical simulations by Cam- 
panelli et al. [s^. We plot the final value of the kick, 
Wkick = v^^'^^{tf), where i/ is the last time in the simula- 
tion, as a function oi /3 — f3o in Fig. [131 The sinusoidal 
dependence in our model is exact up to numerical error. 
One can see this must be the case from examining the 



FIG. 14: The kick as a function of time, for several different 
initial angles of the spin /? — /3o . There is also a straight green 
(light gray) dashed line at zero velocity to indicate how the 
kick varies around this point. 



form of our expression for the momentum flux, Eq. (|103l) . 
Because the evolution equations are not influenced by the 
orientation of the spins, then the Zerilli function will be 
identical for different initial spin directions. The Regge- 
Wheeler function, however, will evolve in the same way, 
but because the value on the matching surface is propor- 
tional to e~*^(*^ [see Eq. ([84])] . the different evolutions 
will also differ by an overall phase, e"^^ , where /3 is the 
initial value of the spin. Thus, when one takes the imag- 
inary part of product of the Regge- Wheeler and Zerilli 
functions to get the momentum flux in Eq. (|103p . one 
will have sinusoidal dependence. (In fact, we could have 
simply done one evolution and changed the phasing as 
described above to find the above results; as a test of 
our method, however, we in fact performed the multiple 
evolutions to confirm this idea.) 

We close this section by making the following obser- 
vation, which may be known, though we have not seen 
in numerical-relativity results. Since the dependence on 
/3 of the kick is sinusoidal, then for each /3, (3 — ir gives 
the same kick and momentum fiux pattern, Pzit), just 
opposite in sign. At the same time, though, because of 
the sinusoidal dependence there are two values that give 
rise to the same kick in the same direction; however the 
shape of the momentum flux Pz{t) is not the same for 
these two. One can see this in Fig. [T31 where the black 
dotted and dashed curve and minus the red (gray) dotted 
and dashed curve give rise to the same kick; nevertheless, 
the pattern of the momentum flux is very different. A 
careful study of this would reveal more about how the 
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FIG. 15: The kick calculated for several initial values of 
/3, minus the angle that produces nearly zero kick, Po — 
2157r/192. The data points are calculated from the numerical 
evolutions of this section, while the solid curve is a sinusoidal 
fit to the data. The hybrid model produces a sinusoidal de- 
pendence of the kick on the initial angle /3 — /3o very precisely. 

spins precess and would be of some interest. 

VI. CONCLUSIONS 

In this paper, we extended a hybrid method for head- 
on mergers to treat inspiralling black-hole binaries. We 
introduced a way to include a radiation-reaction force 
into the hybrid method, and this led to a self-consistent 
set of equations that evolve the reduced-mass motion of 
the binary and its gravitational waves. Using just PN 
and linear BHP theories, we were able to produce a full 
inspiral-merger-ringdown waveform that agrees well in 
phase (though less well in amplitude) with those seen 
in full numerical-relativity simulations. Even though the 
dynamics during inspiral follow the modified dynamics of 
a point particle in Schwarzschild rather than the exact 
dynamics of a black-hole binary, the phasing in the wave- 
form agrees well. Because we assume the background is 
a Schwarzschild black hole (rather than a Kerr, the true 
remnant of black-hole binary inspirals), the merger and 
ringdown parts of the hybrid and numerical-relativity 
waveforms do not match as well. Nevertheless, the hybrid 
method does produce a waveform that is quite similar to 
that of numerical relativity. 

We also studied spinning black holes, particularly the 
superkick configuration (antialigned spins in the orbital 
plane). We discussed a method to incorporate spin 



precession, based on the geodetic precession of a spin- 
ning point particle in the Schwarzschild spacetime. This 
caused the spins to lock to the orbital motion during the 
merger and ringdown, which, in turn, helped to replicate 
the pattern of the momentum flux and the sinusoidal de- 
pendence of the merged black hole's kick velocity seen 
in numerical simulations. Again, because the amplitude 
of the emitted gravitational waves does not match that 
of numerical-relativity simulations, the magnitude of the 
kick does not completely agree. Nevertheless, because 
the approximate method was able to capture the pat- 
tern of the momentum flux, it gives credence to the idea 
the locking of the spin-precession frequency to the orbital 
frequency contributes to large black- hole kicks. 

It would be of interest to extend this approach to see if 
it can recover the results of numerical relativity more pre- 
cisely. To do this would involve a two-pronged approach: 
on the one hand, we would need to include higher PN 
terms in the metric in the interior while using a more 
accurate Hamiltonian to describe the conservative dy- 
namics of the binary (such as the EOB Hamiltonian); on 
the other hand, we would need to evolve the perturba- 
tions in a Kerr background. It would be simpler to choose 
the Kerr background to have the spin of the final, merged 
black hole, but one could also envision evolving perturba- 
tions in an adiabatically changing Kerr-like background 
with a slowly varying mass and angular momentum pa- 
rameter that change in response to the emitted gravita- 
tional waves. It would be of interest to see if such an 
approach leads to an estimate of the spin of the final 
black hole similar to that proposed by Buonanno, Kid- 
der, and Lehner ^6^. Incorporating the PN corrections 
and a new Hamiltonian would be the most straightfor- 
ward improvement, while those involving the Kerr back- 
ground are technically more challenging, and computa- 
tionally more expensive. 
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